clear; close all; clc;


%% 参数设置
alpha = 5;       % 控制参数α
sigma = 0.8;     % 控制参数σ
mu = 0.1;        % 控制参数μ
epsilon = 0.1;   % φ 更新系数
phi0 = 4;        % 忆阻器初始状态

N_total = 2000;  % 总迭代次数
N_cut = 1000;    % 舍去的瞬态迭代次数

% 耦合强度k取值范围
k_min = 0;
k_max = 0.4;
Nk = 1000;
k_values = linspace(k_min, k_max, Nk);

% 用于存储X1和X2的稳态极值
K1 = [];
X1 = [];
K2 = [];
X2 = [];

%% 对每个耦合强度k调用封装函数计算，并提取稳态极值
for idx = 1:length(k_values)
    k = k_values(idx);
    [X1_save, X2_save, E] = Rulkov(alpha, sigma, mu, k, phi0, N_total, epsilon, N_cut);
    
    % 取去除瞬态后的数据
    transient_range = (N_cut+1):N_total;
    X_1 = X1_save(transient_range);
    X_2 = X2_save(transient_range);
    
    % 对 X1 数据提取局部极值（最大值和最小值）
    [peaks1, ~] = findpeaks(X_1);      % 查找局部最大值
    [troughs1, ~] = findpeaks(-X_1);   % 查找局部最小值
    troughs1 = -troughs1;              % 恢复局部最小值
    % 保证行列一致
    peaks1   = peaks1(:);
    troughs1 = troughs1(:);
    extrema1 = sort([peaks1; troughs1]);
    
    % 对 X2 数据提取局部极值
    [peaks2, ~] = findpeaks(X_2);
    [troughs2, ~] = findpeaks(-X_2);
    troughs2 = -troughs2;
    peaks2   = peaks2(:);
    troughs2 = troughs2(:);
    extrema2 = sort([peaks2; troughs2]);
    
    % 当前 k 值对应的数据（分别对 X1、X2）
    K_1 = repmat(k, length(extrema1), 1);
    K_2 = repmat(k, length(extrema2), 1);
    
    % 保存数据
    K1 = [K1; K_1];
    X1 = [X1; extrema1];
    K2 = [K2; K_2];
    X2 = [X2; extrema2];
end

%% 绘制稳态极值分岔图
figure; hold on;

% 绘制 X1 极值 (红色, 有一定透明度)
scatter(K1, X1, 1, 'filled', 'b');

% 绘制 X2 极值 (蓝色, 有一定透明度)
scatter(K2, X2, 1, 'filled', 'r');

xlabel('k');
ylabel('x');
title('分岔图');
legend('X1 极值','X2 极值','Location','best');

hold off;
